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1.0 INTRODUCTION 


The solution of the initial value problem for ordinary differential equations 
(ODE’s) 


dy 

dt 


f(t.y)» 


y(t 0 ) = yo 


(D 


may be treated by several different numerical methods. Runge-Kutta (RK) al- 
gorithms are a type of method well suited to solving equation (1) for many 
classes of functions, f, because of their simplicity and their accuracy. The 
RK algorithm is derived using a direct comparison with a truncated Taylor se- 
ries, giving the accuracy of the Taylor series without the difficulty of deter- 
mining complicated partial derivatives. The comparison between the Taylor se- 
ries expansion of the solution vector and the solution determined by the RK 
algorithm results in a number of expressions referred to as truncation error 
coefficients, T^ j. Associated with each term of order i in the Taylor se- 
ries (or with eacli power of h, the integration stepsize) are Xj. truncation 
error coefficients. For an RK algorithm to be of order p, the Tj^j coeffi- 
cients must be identically zero for i = 1,...,p; j r? l,..,,Xj_, These vanishing 
truncation error coefficients are referred to as equations of condition. The 
nonvanishing error coefficients, however, are of equally great importance since 
they indicate how closely the RK solution approximates a Taylor series solution 
of higher order. The equations of condition determine the validity of an RK 
algorithm; the nonvanishing error coefficients explain the differences between 
particular RK algorithms of the same order. While a user may apply an RK 
algorithm, never considering the truncation error coefficients, awareness of the 
effect of these terms is important both in the selection of a specific algorithm 
and in the analysis of difficulties encountered during the solution of a particu- 
lar ODE. 

D. G. Bettis, a has developed an algorithm for generating the truncation error co- 
efficients for RK methods designed to treat systems of both first- and second- 
order ODE's directly. The recursive nature of this algorithm lends itself 
readily to computer programing, generating high order error coefficients with 
little added difficulty. Such an algorithm, implemented in a numerical code, is 
an essential tool for anyone developing coefficients for RK algorithms and is of 
interest to the user of RK methods in analyzing the effectiveness of specific RK 
algorithms. A Fortran subroutine, RKEQN, written to accompany reference 1, 
generated the truncation error coefficients through order 10 but required a 
great amount of storage location, particularly when a double precision version 
of the program was needed. The basic structure of this original program has 
been reformulated to reduce storage requirements significantly and to accommo- 
date variable dimensioning. This new Fortran program, SUBROUTINE RKEQ, deter- 
mines truncation error coefficients for RK algorithms in the sequence presented 
in reference 1 for orders 1 through 10 and extends the order of coefficients 


a From a private communication with D. G. Bettis, 1978. 
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through 12 with the 11th- and 12th-order terras determined following the patterns 
used to establish the lower order coefficients. Both subroutines (RKEQN and 
RKEQ) are also written to treat RK m-fold methods (refs. 2 and 3) which utilize 
m known derivatives of f to inorease the order of the algorithm. Setting 
m = 0 gives the classical RK algorithm. 
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2.0 THE GENERATION OF TRUNCATION ERROR COEFFICIENTS 

The solution of equation (1) at t-j = t 0 ♦ h ( using the RK algorithm, is written 



( 2 ) 


where 


f 0 - f(t 0 ,y 0 ) 


f k = f(t 0 + a k h, y 0 + h 



&k,X f\) 


where P + 1, the number of evaluations of f computed, is referred to as the 
number of stages. The truncation error coefficients, T-^j, determined by 
comparing the Taylor series expansion of equation (2) wit& the Taylor series ex- 
pansion of the solution about t 0 , are nonlinear combinations of the C, a, and 
8 coefficients. For the classical RK algorithm, the jth error coefficient of 
order i assumes the form 

T l,j = {*1,J - %,J (3) 

j = 1,...,Xj_ where p = i. For the m-fold algorithm, p = m -f 1 and corre- 
sponds to the order of the term. (For m = 0, the m-fold algorithm is identical 
to the classical RK formula.) The A^j and j terms are constants (or 
functions of m for m-fold methods) ana may be determined by recursive 
relations. (One should note that while references to m-fold RK algorithms may 
appear to complicate matters, the inclusion of these methods in RKEQ (or RKEQN) 
involves the insertion of only a few additional lines of coding. Once these ad- 
ditions are made, the classical RK error coefficients and the m-fold error coef- 
ficients are determined identically.) The complicated expression to generate in 
equation (3) is the F-j^j term, which is a combination of the C, a, and 8 
coefficients ’ 


P 

F i,j = S ° k Si >3> k ‘ ^ 

k=k 0 

with k being a combination of a and 8 coefficients and where k 0 de- 

pends upon’ the number of summations embedded in j k . 
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The algorithm developed by Bettis and us«)d for generating these S^j^, A^j, 
and B^ j terms is documented in reference 1 , where the u terms are 

writtenin an abbreviated notation, e.g., Sq « = a 2 0a00a, with the subscript- 
ing and embedded summations being suppressed. The rules for writing the entire 
Si j terms are also described thoroughly in reference 1. For the sake of 
Interpreting the program, however, a few features need to be known about gener- 
ating the terms in abbreviated notation i Denoting the number of truncation 
error coefficients of order p, by Ay, and suppressing the k subscript, the 
first 2A u -i S^ j u expressions are H generated from the fcerms * The 

remaining H An - expressions, referred to as composite sums, are formu- 
lated as products of lower order S terms. The and constants are 

also generated from simple relationships involving previous A and B terms. 

In generating the first Aj,_i terms of order i, the expressions are 

preraultlplied by an a. (Adjacent a's represent actual multiplication.) Thus, 
Sg 13 = aSQ, S 13 = a 3 0a00ot. The next Ai_i terms, S if j, j s Xi_i + 

1 ,..., 2 X^_i are generated by premultiplying the & i- 1 , j-Aj-i expressions by a 
3, e.g., Sg i 20 = BSq 13 = 0a 2 0a00a • (Adjaoent 0's do not represent multipli- 
cation of the 0 coefficients since a summation sign precedes each 0 when the 
entire S term is written, e.g., 


s 6,4,k = <* 2 33ot = <*k 2 




3ji,rrf*m 


(5) 


( s 6,4,k * ot 2 0 2 ot.)) The recursion relations, then, for the first Ai_i terms 
of order 1 , (i > 2 ), are 

s i,j s aS i- 1 , j 
A l,j = 

B i,j = M B i- 1 ,j 


for j = 1 , 2 , . . . , X i_ i , where p is the power of the leading a in the 
term, and 


s i,j = 3 s i-1»J-A i _ 1 


Ai,j = (m+i- 1 ) A i . 1>j _ x ._ i 


4 
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B i»j s 


for j s + 1 , . . . » 2 Xi_i • (Sij = 1 , and S2, 1 « a.) The first Xi»i S 
expressions are referred to as alpha terms in subroutine RKEQ, while the next 
Xi_i expressions are called beta terms. The remaining terms, composite sums, 
are generated by considering the weight factors of the S terms. The Sp ( j ex- 
pressions have a weight factor of y-1 (l.e., the number of a and 3 coeffi- 

cients Included in the S term) . The composite sums of order i are all prod- 
ucts of S^j terms having initial 3 coefficients, whose weight factors add 
up to i- 1 . Subroutine RKEQ determines these composite sums in a separate block 
of the subprogram, calling subroutine CROSS to perform the multiplication of the 
S, A, and B terms. The j and B^j terms of a composite sum are the 
products of the A and B constants whose corresponding S terms form the com- 
posite sum. (When an S term is raised to the power k, an additional kl 
multiplies the B constant.) 


5 
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3.° DESCRIPTION OF THE FORTRAN PROGRAM 

Subroutine RKEQ determines the truncation error coefficients (TEC) for a given 
set of RK coefficients and returns TERROR, 

9 )l/2 

TERROR * 2, T i,J 


for a specified order i. Since RK algorithms with embedded pairs of solutions, 
e.g. RK-Fehlberg formulas, are often studied, RKEQ is written to treat two 
algorithms simultaneously, which use identical a and $ coefficients but 

A A 

different and C k coefficients. (TERROR is formed using the C^ 
coefficients.) The Greek letters and 8^ are replaced by A(I) and B(I) and 
the Aj^j and B-^j constants are denoted AA(I,J) and BB(I,J), respectively. 

The input parameters for RKEQ are; 

(1) The RK coefficients 

(a) A(K) cty, the alpha coefficients 

(b) CO, C(K) C 0 , Cfc, the C^ coefficients for the first solution 

A A A 

(c) CHO, CH(K) C 0 , C^, the C^ coefficients for the second solution 

used to form TERROR 

(d) BO(K) , B(K,L) 8k, o> 3k the beta coefficients 


where K = 1,2,...,R, L = 1,2,..., K - 1, R an integer with R + 1 being the 
number of stages of the algorithm, and 

(2) The integers controlling orders and options 

(a) R = the index for dimensioning the RK coefficients 

(b) IORDER = the maximum order to be treated 

(c) ITERR = the order of TEC used to form TERROR 

(d) IOPT = the options for operating the program. For IOPT = 1, RKEQ com 
putes and prints all TEC(I,J) for I s 1 , . . . , IORDER. For IOPT = 2, 
RKEQ computes and prints TEC( IORDER, J) only. For IOPT = 3, RKEQ com 
putes but does not print TEC( ITERR, J). (For all options TERROR is 
computed, which may require internal adjustments to the order.) 

(e) MFOLD = an integer giving the number of known derivatives of f. For 
the classical RK algorithm, MFOLD = 0. 
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(f) LS a a dimensioning index for the work arrays, S, AA and BB, 

The output parameter for RKEQ is TERROR, the Euclidean norm of the TEC of order 
ITERR, Depending upon the option used, RKEQ may print values of TEC, but these 
are not returned to the main program. 

Parameters S, AA, and BB are used internally by RK'vj to compute the TEC 
terms. To take advantage of variable dimensioning, these parameters are given 
in the calling sequence with dimensions S(LS,R), AA(LS), BB(LS) . The RK coef- 
ficients should be dimensioned A(R) , C(R), CH(R), B(R, R), B0(R), R an inte- 
ger, where R + 1 is the number of stages for the algorithm, 

The calling sequence for RKEQ is 

SUBROUTINE RKEQ(A, C, CO, CH, CHO, B, BO, R, IORDER, ITERR, IOPT, MFOLD, 
TERROR, S, LS, AA, BB). 

which, if a printing option is used, will give the TEC from the C solution 
in the first column and the TEC from the CH solution in the second column. 
Integers IORDER, ITERR, and IOPF <?re reset within the subroutine, and any adjust- 
ments made to protect against exceeding dimension or option limits are made to 
the new variables, so that the user may enter constant values in the calling 
sequence of the driving program. 

A sample calling sequenoe for a six-stage, fifth-order algorithm is 

CALL RKEQ(A, C, CO, CH, CHO, B, BO, 5, 7, 6, 1, 0, TERROR, S, 48, AA, BB) 

which computes and prints all TEC through order 7 for the classical RK formulas, 
using the 6th-order terms of the CH solution to form TERROR. Using the calling 
sequence 


CALL RKEQ(A, CH, CHO, C, CO, B, BO, 5, 7, 6, 1, 0, TERROR, S, 48, AA, BB) 
generates similar information exoept that TERROR is formed by the C solution 

A 

(and the C TEC terras are now printed in the first column.) 

The minimum value of LS for a given ORDER, I, is found in table I. A listing 
of subroutines RKEQ and CROSS may be found in the appendix. 


TABLE I.- DIMENSIONING PARAMETER, LS, FOR ORDERS 1 THROUGH 12 


ORDER 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

LS M in 

1 

1 

2 

A 

9 

20 

48 

115 

286 

719 

1842 

4766 


7 
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H.O CONCLUDING REMARKS 


A program, which evaluates thi* truncation error coefficients, is an essential 
tool in the development of Runge-Kutta algorithms and in the comparison of 
existing RK algorithms. By structuring the routine in the given form, a substan- 
tial savings in at ^rage occurs in generating these truncation error coefficients 
using the recursive formulation presented by D. G. Bettis (ref. 1). The exten- 
sion to orders higher than 12 is relatively simple but not of great practical 
use at the present time. 
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APPENDIX A 


SUBROUTINE RKEQ 


S0FM1 3 


SUBROUTINE RKEQ ( A , C , CO , CH , CHO , B , BO , R , 

1 I ORDER , ITERR , IOPT ,MFOLD .TERROR ,S , LS , AA ,BB) 

INTEGER R, ORDER, OPTION 

DOUBLE PRECISION A(R) ,C(R) ,CH(R) ,B(R,R) ,BO(R) 

DOUBLE PRECISION CO CHO 

DOUBLE PRECISION S(LS,R) ,AA(LS) ,BB(LS) 

DOUBLE PRECISION SUM1 ,SUM3, XI , ZERO, UNITY, TWO, 

1 MP(12) ,MM(12) , FACT (12) 

DOUBLE PRECISION P ,PP1 ,PM1 ,PM2 ,ZAPP .TERROR 
DIMENSION M( 12) , LIMIT (12), INDEX (11) 

LOGICAL EVEN 
C 

DATA LIMIT/1,1,2,4,9,20,48,115,286,719,1842,4766/ 

DATA INDEX/1,1,2,5,11 ,28,67,171 ,433,1123,2924/ 

DATA ZERO , UNITY , TWO , ZAPP/O . ODO , 1 . ODO , 2 . ODO , 1 . OD- 1 4 / 
DATA MAXORD/12/ 

C 

c 

DATA KPRINT/O/ 

C 

QiOOftttftftttftttttftttttttftfttfftttftttttilttftftft 

C 

C 

C 

C SUBROUTINE RKEQ IS A FORTRAN SUBROUTINE WRITTEN BY 

C M.K. HORN WHICH IMPLEMENTS THE ALGORITHM DEVELOPED 

C BY D.G. BETTIS TO GENERATE TRUNCATION ERROR COEFFI- 
C IENTS FOR RUNGE-KUTTA ALGORITHMS. 

C 

C REFERENCE: BETTIS, D.G. AND M.K. HORN, * COMPUTATION 
C OF TRUNCATION ERROR TERMS FOR RUNGE-KUTTA METHODS,’ 

C TICOM REPORT 77-14, DECEMBER, 1978. 

C 

C SUBROUTINE RKEQ DETERMINES THE TRUNCATION ERROR 
C COEFFICIENTS (TEC) FOR A RUNGE-KUTTA ALGORITHM HAVING 
C AN EMBEDDED PAIR OF SOLUTIONS 

C 

C R 

C 

C Y = Y + H SUM C F 
CIO K K 

C KsO 

C 
C 

C R 

C 

C YH = Y + H SUM CH F 
C 10 K K 


00000200 

00000300 

00000400 

00000500 

00000600 

00000700 

00000800 

00000900 

00001000 

00001100 

00001200 

00001300 

00001400 

00001500 

00001600 

00001700 

00001800 

00001900 

00002000 

0000210O 

00002200 

00002300 

00002400 

00002500 

00002600 

00002700 

00002800 

00002900 

00003000 

00003100 

00003200 

00003300 

00003400 

00003500 

00003600 

00003700 

00003800 

00003900 

00004000 

00004100 

00004200 

00004300 

00004400 

00004500 

00004600 

00004700 

00004800 

00004900 


ooooooooooooooooooooooooooooooooooo 


00005000 


C 


K=0 


WHERE F = F(T ,Y ) 
0 0 0 


K-1 


F = F(T + A H, Y + H » SUM B F ) 
K 0 K 0 K,L L 

L=0 


FOR MFOLD = 0, RKEQ FORMS THE TEC FOR THE CLASSICAL RK 
FORMULAS OF ORDER = IORDER. FOR MFOLD = M, RKEQ FORMS 
THE TEC FOR MFOLD-RK FORMULAS OF ORDER = M + IORDER. 


OPTION .EQ. 1 

OPTION .EQ. 2 
OPTION .EQ. 3 


COMPUTES AND PRINTS ALL TRUNCATION 
ERROR COEFFICIENTS THROUGH ORDER = 
IORDER 

COMPUTES AND PRINTS ONLY T.E.C. OF 
ORDER = IORDER 

COMPUTES T.E.C. AS IN OPTION = 2 BUT 
DOES NOT PRINT 


ADJUSTS INPUT PARAMETERS IF THESE ARE NOT WITHIN 
ALLOWABLE RANGE 


ORDER = IORDER 
LIMITS = ITERR 

IF (ORDER .LT. LIMITS) ORDER = LIMITS 
IF (LS. LT. LIMIT (MAXORD)) GO TO 1 
IF (ORDER .LE. MAXORD) GO TO 1 
ORDER = MAXORD 
PRINT 507, MAXORD 

507 FORMAT (52H ORDER REQUESTED IS BEYOND CAPABILITY OF THE PROGRAM 
1 ,/,24H ORDER LOWERED— ORDER = ,12) 

IF (ITERR .LE. ORDER) GO TO 1 
LIMITS = ORDER 
PRINT 526, LIMITS 

526 FORMAT (45H ITERR REQUESTED IS LARGER THAN MAXIMUM ORDER 
1 ,/,35H ITERR HAS BEEN REDUCED TO ORDER = ,12 ) 


00005100 

00005200 

00005300 

00005400 

00005500 

00005600 

00005700 

00005800 

00005900 

00006000 

00006100 

00006200 

00006300 

00006400 

00006500 

00006600 

00006700 

00006800 

00006900 

00007000 

00007100 

00007200 

00007300 

00007400 

00007500 

00007600 

00007700 

00007800 

00007900 

00008000 

00008100 

00008200 

00008300 

00008400 

00008500 

00008600 

00008700 

00008800 

00008900 

00009000 

00009100 

00009200 

00009300 

00009400 

00009500 

00009600 

00009700 

00009800 

00009900 


o o o 


C 


00010000 


1 CONTINUE 
C 

TERROR = ZERO 
OPTION = IOPT 

IF (IOPT , LE. 0 .OR. IOPT ,GT. 3) GPT.U.N = 1 
IF (LS .GE. LIMIT (ORDER)) GO TO 3 

2 CONTINUE 

ORDER = ORDER - 1 

IF (LS .LT. LIMIT (ORDER)) GO TO 2 
PRINT 505 , ORDER , IORDER , LIMIT ( IORDER ) 

505 FORMAT (50H THE ORDER SPECIFIED HAS BEEN REDUCED TO ORDER = 

1 ,I3,/,32H BECAUSE OF INSUFFICIENT STORAGE ,/, 

2 14H FOR ORDER = ,I3,20H LS MUST BE .GE. ,15, 

3 /,56H TERROR s SQRT(SUM(T.E.C. (I,J)*T.E.C.(I,J))) IS COMPUTED 

4 ,/, 14H FOR I = ORDER ) 

C 

IF (LIMITS .LE. ORDER) GO TO 3 
LIMITS = ORDER 
PRINT 527, LIMITS 

527 FORMAT (45H ITERR REQUESTED IS LARGER THAN THE PROVIDED 
1 ,/,50H DIMENSIONING. ITERR HAS BEEN REDUCED TO ITERR = , 

2 12) 

3 CONTINUE 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 


IF (ORDER .GT. MAXORD) ORDER = MAXORD 


MP( J) = (MFOLD+J) D.P. 

M(J) = J INTEGER 

FACT(J) = J D.P. 

MM( J) = (MFOLD+J) D.P. 


INDEX— COUNTS THE NUMBER OF S(J,K) WITH A GIVEN A«*K 
AS THE FIRST TERM IN THE EXPRESSION, E.G., 

FOR J=7, THERE ARE 1 A**6, 1 A*#4, 2 A**3, 

4 A**2, AND 9 A»*1 §S 


XI = A( 1 ) - BO ( 1 ) 

JJ = 1 

550 FORMAT ( 1 5H ERROR IN BETA( ,12, 10H) SUM = ,D15.7) 

IF (DABS (XI ) .GE. ZAPP) PRINT 550,JJ,X1 
DO 5 J = 2,R 
XI = A( J) - B0( J) 


00010100 

00010200 

00010300 

00010400 

00010500 

00010600 

00010700 

00010800 

00010900 

00011000 

00011100 

00011200 

00011300 

00011400 

00011500 

00011600 

00011700 

00011800 

00011900 

00012000 

00012100 

00012200 

00012300 

00012400 

00012500 

00012600 

00012700 

00012800 

00012900 

00013000 

00013100 

00013200 

00013300 

00013400 

00013500 

00013600 

00013700 

00013800 

00013900 

00014000 

00014100 

00014200 

00014300 

00014400 

00014500 

00014600 

00014700 

00014800 

00014900 


A- 5 


amt 




o o o o o o oooooooo 


JLOW = J-1 
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00015000 


DO U : 1 , JLOW 

4 XI s XI - B( J ,K) 

IF (DABS (XI ) .GE. ZAPP) PRINT 550, J, XI 

5 CONTINUE 
C 

T.F (KPRINT .EQ. 1) READ PRINT 499 


SETS VALUES OF FACTORIALS USED IN THE PROGRAM. 
ADJUSTMENTS IN THE PROGRAM TO ACCOMODATE M-FOLD 
RUNGE-KUTTA ALGORITHMS OCCUR HERE. IF MFOLD = 0, 
THE CLASSICAL RUNGE-KUTTA T.E.C. OCCUR 


IF (MFOLD .GT. 0) GO TO 6 
EVEN = .TRUE. 

GO TO 8 
6 CONTINUE 

N1 = MFOLD / 2 
N2 a 2 • N1 
EVEN = .FALSE. 

IF (N2 .EQ. MFOLD) EVEN = .TRUE. 
8 CONTINUE 

MP( 1 ) = DFLOAT( MFOLD )+UNITY 
M(1 ) = 1 
XI = UNITY , 

FACT(1 ) = UNITY 
DO 10 I = 2, ORDER 
XI s XI + UNITY 
FACT(I) = FACT(I-1)»X1 
MP( I) = MP(I-I) + UNITY 

10 M(I) = M(I-1) + 1 
M2 = 1 

Ml = 0 

11 CONTINUE 

Ml = Ml + 1 
M2 s M2*M1 

IF (Ml .LT. MFOLD+1 ) GO TO 11 
MM( 1 ) = DFLOAT(M2) 

DO 12 I = 2, ORDER 

12 MM(I) = MP(I)» MM(I-I) 


IF (KPRINT .EQ. 1) PRINT 499 

SETS INITIAL VALUES OF S(I,J) EQUAL TO ZERO 
AND SETS AA( J) AND BB(J) EQUAL TO UNITY 

DO 14 J = 1 , LS 


00015100 

00015200 

00015300 

00015400 

00015500 

00015600 

00015700 

00015800 

00015900 

00016000 

00016100 

00016200 

00016300 

00016400 

00016500 

00016600 

00016700 

00016800 

00016900 

00017000 

00017100 

00017200 

00017300 

00017400 

00017500 

00017600 

00017700 

00017800 

00017900 

00018000 

00018100 

00018200 

00018300 

00018400 

00018500 

00018600 

00018700 

00018800 

00018900 

00019000 

00019100 

00019200 

00019300 

00019400 

00019500 

00019600 

00019700 

00019800 

00019900 


A-6 


OOOO OOOOO O O o o 


AA ( J ) = UNITY 


00020000 


BB( J) = UNITY 


00020100 

DO 14 K a 1,R 


00020200 

14 S(J,K) a ZERO 


00020300 

00020400 

IF (OPTION .GT. 1 .AND. LIMITS ,NE. 1) GO TO 21 

00020500 

00020600 

EVALUATES T.E.C, OF ORDER 1 


00020700 

00020800 

SUM1 = UNITY/DFLOAT (MFOLD + 

1) - CO 

00020900 

SUM3 = UN ITY /DFLOAT ( MFOLD + 

1) - CHO 

00021000 

DO 20 I = 1,R 


00021100 

SUM1 = SUM1 - C(I) 


00021200 

SUM3 = SUM3 - CH(I) 


00021300 

20 CONTINUE 


00021400 

00021500 

JJ1 = 1 


00021600 

IF (LIMITS ,EQ. 1) TERROR a 

SUM3*SUM3 

00021700 

IF (OPTION .EQ. 3) GO TO 21 


00021800 

PRINT 500, JJ1 


00021900 

PRINT 501 , JJ1 , JJ1 ,SUM1 ,SUM3 
500 FORMAT (36H TRUNCATION ERROR 


00022000 

TERMS X-ORDER s ,I2,2X 

00022100 

1 ,//) 


00022200 

501 FORMAT (2 (2X,I4) ,2(D15.7)) 


00022300 

00022400 

21 CONTINUE 


00022500 

00022600 

SETS S(1,J) TERMS 


00022700 

00022800 

KOUNT = 1 


00022900 

00023000 

DO 26 I = 1,R 


00023100 

IF (MFOLD ,GT. 0) GO TO 22 


00023200 

S(1 ,1) = A(I) 


00023300 

GO TO 26 


00023400 

22 CONTINUE 


00023500 

IF (EVEN) GO TO 24 


00023600 

S(1,I) = DABS (A (I) )##MP( 1 ) 


00023700 

GO TO 26 


00023800 

24 CONTINUE 


00023900 

XI a DABS ( A ( I ) ) ##MP ( 1 ) 


00024000 

S(1 ,1) = DSIGN(X1,A(I)) 


00024100 

26 CONTINUE 


00024200 

00024300 

EVALUATES T.E.C. FOR ORDER = 

KOUNT =2,3 

00024400 

00024500 

00024600 

AA(1) = UNITY 


00024700 

BB( 1 ) = MM( 1 ) 


00024800 

28 CONTINUE 

A- 7 

00024900 
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P S TWO + DFLOAT(MFOLD) 


00025000 


PP1 = P+UNITY 
30 CONTINUE 

IF (LIMITS ,EQ. KOUNT) JJ1 = K0UNT - 1 

IF (OPTION ,GT, 1 .AND. LIMITS .NE. KOUNT) GO TO 34 

JJ2 = 0 

JJ1 = KOUNT + 1 

PRINT 500, JJ1 

DO 33 K s 1, KOUNT 

JJ2 s JJ2 + 1 

SUM1 = UNITY/ (AA(K)«P) 

SUM3 = SUM1 
DO 32 I = 1,R 

SUM1 = SUM1 - C(I)*S(K,I) 

SUM3 = SUM3 - CH(I)*S(K,I) 

32 CONTINUE 

IF (LIMITS .EQ. KOUNT) TERROR = TERROR + SUM3*SUM3 
IF (OPTION .EQ. 3) GO TO 33 
PRINT 501 , JJ1 , JJ2,SUM1 ,SUM3 

33 CONTINUE 

34 CONTINUE 

IF (KOUNT .EQ. 2) GO TO 38 

SETS S(1,J), S(2,J) FOR THIRD ORDER T.E.C. 

KOUNT = 2 
P = P + UNITY 
PP1 = PP1 + UNITY 
DO 35 I = 2,R 
S(2 , I) = ZERO 
IM1 = 1-1 
DO 35 J = 1 t IM1 

35 S(2,I) = S(2,I) + B(I , J)*S( 1 ,J) 

DO 36 I = 1,R 

36 S( 1 ,1) = S(1,I)»A(I) 

AA(2) = TWO 

BB(1 ) = MM(2) 

BB(2) = MM( 1 ) 

GO TO 30 
38 CONTINUE 

IF (ORDER .LE. 3) GO TO 182 
LIM1 = 3 

EVALUATES T.E.C. FOR ORDERS GREATER THAN THREE 

DO 180 J = 4, ORDER 
P = P + UNITY 
PP1 = P + UNITY 
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PM1 s P - UNITY 


00030000 


PM2 s P - TWO 

IF (KPRINT ,EQ. 1) READ 1497, II 
LIM1 s LIM1 + 1 
LIMA s LIMIT(J-I) 

LIMB s LIMA 

COMPUTES S(1,I)— ALL OTHER S(J,K) TERMS INVOLVING A 

LEADING ALPHA ARE ALREADY DETERMINED 
EXCEPT FOR THE POWER OF ALPHA WHICH 
IS DETERMINED BY IN INDEX AND J 


AA ( 1 ) s AA( 1 ) 

BB( 1 ) s BB(1) • MP(LIMI-I) 
DO 42 I s 1,R 
42 S<1 f I) s S(1,I)»A(I) 

LIM1 s INTEGER P 


BETA TERMS 

MARKB = LIMA+1 
DO 61 I s 2,R 
S( MARKB, I) = ZERO 
IM1 s 1-1 
DO 61 K = 1 , IM1 

61 S(MARKB,I) = S ( MARKB , I ) +B ( I , K ) * A ( K ) ** ( MFOLD+LIM1 -2 ) 
AA (MARKB )sPM1 
BB( MARKB) = MM(LIM1-2) 

IND2 = 1 

IND1 s MARKB 

LL = LIM1 - 4 

IF (LIM1 .EQ. 4) GO TO 66 

DO 65 K r. 1,LL 

LL1 = INDEXU+1) 

IPOW = LL-K+1 
DO 65 KK = I ,LL1 
IND1 = IND1 + 1 
IND2 = IND2 + 1 
DO 64 I = 2,R 
S(IND1,I) = ZERO 
IM1 = 1-1 
DO 64 L = 1 , IM1 

64 S(IND1,I) a S(IND1,I)+B(I,L)»S(IND2,L)»A(L)«M(IP0W) 
BB(INDI) = BB ( IND2 ) * FACT (IPOW ) 

65 AA(INDI) = AA(IND2)#PM1 

66 CONTINUE 

LL 8 LIMB - IND2 
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DO 68 K = 1,LL 


00035000 


IND1 s IND1 + 1 
IND2 s IND2 + 1 
DO 67 I = 2,R 
S(IND1,I) n ZERO 
IM1 = 1-1 
DO 67 L a 1 ,IM1 

67 S(IND1 ,1) = S(IND1,I)+B(I,L)*S(IND2,L) 

BB(INDI) = BB(IND2) 

68 AA<IND1) = AA(IND2)*PM1 
883 CONTINUE 

JM3 = J - 3 

GO TO (150, 100, 105, 110, 115, 120, 125, 130, 135), JM3 


CROSS PRODUCT TERMS 
100 CONTINUE 

5 TH ORDER TERMS 
IND = 2«LIMIT(4)+1 

CALL CR0SS(LS,R,S,AA,BB,IND,2,2,0,0) 
GO TO 150 

6 TH ORDER TERMS 
105 CONTINUE 

IND = 2*LIMIT(5)+1 

CALL CROSS(LS,R,S,AA,BB,IND,2,1,3,4) 
GO TO 150 


7 TH ORDER TERMS 

110 CONTINUE 

IND = 2*LIMIT(6)+1 

CALL CR0SS(LS,R,S,AA,BB,IND,2,1,5,8) 
CALL CR0SS(LS,R,S,AA,BB,IND,2,3,0,0) 
CALL CR0SS(LS,R,S,AA,BB,IND,3,2,0,0) 
CALL CR0SS(LS,R,S,AA,BB,IND,4, 2,0,0) 
CALL CROSS(LS,R,S,AA,BB,IND,3, 1,4,4) 

GO TO 150 

8 TH ORDER TERMS 

115 CONTINUE 

IND = 2*LIMIT(7 )+1 

CALL CROSS(LS,R,S, AA,BB,IND,2 ,1,10,18) 
CALL CROSS ( LS , R , S , AA , BB , IND , 3'» 1 , 5 , 8 ) 
CALL CROSS(LS,R,S,AA,BB,IND,4, 1,5,8) 
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00040000 


CALL CROSS ( LS , R , S , AA , BB , IND , 2 , 2 , 3 , 4 ) 


GO TO 150 
C 

C 9 TH ORDER TERMS 

120 CONTINUE 
IND s 2*LIMIT (8)4-1 

CALL CR0SS(LS,R,S,AA,BB,IND,2,1,21,40) 
CALL CROSS(LS,R,S,AA,BB, IND, 3, 1,10,10) 
CALL CROSS(LS,R,S, AA,BB,IND,4, 1,10,18) 

DO 121 K = 5,8 

121 CALL CROSS(LS,R,S,AA,BB,IND,K,1 ,K,8) 

CALL CROSS(LS,R,S,AA,BB,IND,9,1,5,8) 

CALL CROSS ( LS , R , S , AA , BB , IND , 2 , 4 , 0 , 0 ) 

CALL CROSS(LS,R,S,AA,BB,IND,2,1,46,48) 

GO TO 150 

10 TK ORDER TERMS 

125 CONTINUE 
IND s 2»LIMIT(9) 4- 1 

CALL CROSS (LS , R, S , AA , BB , IND , 2 , 1 , 49 , 96 ) 
CALL CROSS ( LS , R , 3 , AA , BB , IND , 3 , 1 , 2 1 , 40 ) 
CALL CROSS(LS,R,S,AA,BB,IND,4,1,21,40) 

DO 126 K = 5,9 

126 CALL CROSS(LS, R,S,AA,BB, IND, K, 1,10,18) 
CALL CROSS ( LS , R , S , AA , BB , IND ,2,1,106,113) 
CALL CROSS ( LS , R , S , AA , BB , IND , 2 , 3 , 3 , 4 ) 

CALL CROSS ( LS , R , S , AA , BB , IND , 3 , 3 , 0 , 0 ) 

CALL CR05S(LS,R,S,AA,BB,IND,4,3,0,0) 

CALL CROSS ( LS , R , S , AA , BB , IND , 4 , 2 , 3 , 3 ) 

CALL CROSS(LS,R,S,AA,BB,IND,3,2,4,4) 

GO TO 150 

130 CONTINUE 

11TH ORDER TERMS 

IND s 2*L1MIT(10) 4- 1 

CALL CROSS(LS,R,S,AA,BB,IND,2,1,116,230) 
CALL CROSS ( LS , R , S , AA , BB , IND , 3 , 1 ,49,96) 
CALL CROSS(LS,R,S,AA,BB,IND,4,1,49,96) 

DO 131 K b 5,8 

131 CALL CROSS (LS , R , S , AA , BB , IND , K , 1 , 2 1 , 40 ) 

DO 132 K = 10,18 

132 CALL CROSS ( LS , R , S , AA , BB , IND , K , 1 , K , 1 8 ) 
CALL CROSS ( LS , R , S , AA , BB , IND ,2,2,21,40) 
CALL CR0SS(LS,R,S,AA,BB,IND,19,1,10,18) 
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00045000 


CALL CROSS (LS, R,S,AA,BB, IND, 20 ,1,10,18) 


CALL CROSS(LS,R,S,AA,BB, IND, 2, 1,269,278) 
DO 133 K = 46,48 

133 CALL CROSS(LS,R,S,AA,BB,IND,K,1,5,8) 

CALL CROSS(LS,R,S,AA,BB,IND,2,3,5,8) 

CALL CROSS ( LS , R ,S , AA , BB , IND ,2,2,46,48) 
CALL CROSS ( LS , R , S , AA , BB , IND , 2 , 5 , 0 , 0 ) 

GO TO 150 

12TH ORDER TERMS 

135 CONTINUE 

IND s 2*LIMIT( 1 1 ) + 1 


CALL CROSS(LS,R,S,AA,BB,IND,2,1,287,572) 

CALL CROSS(LS,R,S,AA,BB,IND,3,1,116,230) 

CALL CROSS ( LS , R , S , AA , BB , IND , 4 , 1 , 1 1 6 , 230 ) 

DO 141 K = 5,8 

141 CALL CROSS (LS,R,S,AA,BB, IND, K, 1,49,96) 

DO 142 K = 10,18 

142 CALL CROSS(LS,R,S,AA,BB,IND,K, 1 ,21 ,40) 

CALL CROSS ( LS , R , S , AA , BB , IND ,2,2,49,96) 

CALL CROSS(LS,R,S,AA,BB,IND, 19,1,21 ,40) 

CALL CROSS(LS,R,S,AA,BB, IND, 20,1 ,21,40) 

DO 143 K = 41,48 

143 CALL CROSS(LS,R,S,AA,BB, IND,K, 1 ,10,18) 

CALL CR0SS(LS,R,S,AA,BB,IND,114,1,5,8) 

CALL CROSS(LS,R,S,AA,BB,IND,115,1,5,8) 

DO 144 K = 269,278 

144 CALL CROSS(LS,R,S,AA,BB,IND,K,1,3,4) 

CALL CROSS (LS , R , S , AA « BB , IND , 3 , 3 , 2 , 2 ) 

CALL CROSS(LS,R,S, AA,BB,IND,4,3|2,2) 

CALL CROSS ( LS , R , S , AA , BB , IND ,3,2,20,20) 

CALL CROSS(LS,R,S,AA,BB, IND,4 ,2, 19,19) 

CALL CROSS(LS,R,S,AA,BB,IND,2,4,3,4) 

150 CONTINUE 

TEMPORARY INSERT TO CHECK VALUES OF AA AND BB COEFF 

LL = LIMIT (LIM1) 

IFAKE = 0 

DO 153 K = 1 ,LL 

IFAKE = IFAKE + 1 

IF (IFAKE ,LT. 40) GO TO 153 

IF (KPRINT .EQ. 1) READ 497, II 

IFAKE = 0 

153 PRINT 506,K,AA(K),K,BB(K) 
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C 1 I4,4H) C ,D15.7) 

IF (KPRINT ,EQ. 1) PRINT 499 
C 
C 
C 

IF (LIMITS ,EQ. J) JJ1 s J - 1 
IF (OPTION .GT. 1 ,AND. LIMITS .NE. J) GO TO 100 
C 

JJ2 a 1 
JJ1 E J 
C 

C EVALUATES FIRST T.E.C. OF ORDER J 
C 

SUM1 a UNITY/(AA(1)«P) 

SUM3 = SUM1 

DO 154 I a 1,R 

SUM1 a SUM1 - C(I)«S(1,I) 

SUM3 a SUM3 - CH(I)*S(1,I) 

154 CONTINUE 

SUM1 a SUM1 / BB( 1 ) 

SUM3 a SUM3 / BB( 1 ) 

IF (LIMITS .EQ. J) TERROR s TERROR + 8UM3*SUM3 
IF (OPTION .EQ. 3) GO TO 155 
PRINT 501,JJ1,JJ2,SUM1,SUM3 

155 CONTINUE 
C 

C 

C EVALUATES T.E.C. FOR S(I,J) TERMS WITH ALPHA AS 
C LEADING COEFFICIENT (I .NE. 2) 

C 

IFAKE a 1 
K a 1 
KNT a 2 

IPOW a LIM1 - 3 
LIMD a INDEX (KNT) 

156 CONTINUE 

DO 159 KK = 1.LIMD 
K a K + 1 

SUM1 a UNITY/(AA(K)*P) 

SUM3 = SUM1 

DO 157 I a 1,R 

SUM1 = SUM1 - C(I)»A(I)WlPOW»S(K f I) 

SUM3 = SUM3 - CH(I)«A(I)**IPOW»S(K t I) 

157 CONTINUE 

SUM1 a SUM1 / (BB(K)*FACT(IPOW)) 

SUM3 a SUM3 / (BB(K)*FACT(IPOW) ) 

IF (LIMITS .EQ. J) TERROR a TERROR SUM3*SUM3 
IF (OPTION .EQ. 3) GO TO 159 ' 

JJ2 a JJ2 + 1 
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IF (IFAKE .LT. *10) 00 TO 158 

IFAKE s 0 

IF (KPRINT .EQ. 1) PRINT 499 

158 CONTINUE 
497 FORMAT (13) 

PRINT 501,JJ1,JJ2,SUM1,SUM3 

159 CONTINUE 
C 

KNT = KNT + 1 
LIMD s INDEX(KNT) 

IPOW s IPOW - 1 
IF (IPOW .GE. 1) GO TO 156 
C 

LIMD s LIMIT (LIM1) - LIMA 
C 

C EVALUATES T.E.C. FOR S(I,J) TERMS WITH BETA AS 
C LEADING COEFFICIENT 
C 

DO 162 KK s 1 ,LIMD 
K s K + 1 

SUM1 s UNITY/(AA(K)*P) 

SUM3 s SUM1 

DO 160 I s 1 ,R 

SUM1 = SUM1 - C(I)*S(K,I) 

SUM3 = SUM 3 - CH(I)»S(K,I) 

160 CONTINUE 

SUM1 = SUM1 / BB(K) 

SUM3 = SUM 3 / BB(K) 

IF (LIMITS .EQ. J) TERROR s TERROR + SUM3*SUM3 

IF (OPTION .EQ. 3) GO TO 162 

IFAKE = IFAKE + 1 

IF (IFAKE .LT. 40) GO TO 161 

IFAKE = 0 

IF (KPRINT .EQ, 1) PRINT 499 

161 CONTINUE 

JJ2 = JJ2 + 1 

C PRINT 556,K,LIM1 ,AA(K) ,P,PP1 
556 FORMAT ( 2 (2X ,14 ),3(2X,D15.7)) 

PRINT 501 ,JJ1 ,JJ2,SUM1 ,SUM3 

162 CONTINUE 

IF (KPRINT .EQ. 1) PRINT 499 
499 FORMAT (//) 

180 CONTINUE 
182 CONTINUE 

TERROR = DSQRT( TERROR) 

RETURN 

END 
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SUBROUTINE CR0SS( LS , R ,S , AA ,BB , INDEX ,TERM1 } POWER , TERMS , 
1 TERM3) 

INTEGER TERM1, TERMS, TERM3, POWER, R 

DOUBLE PRECISION S(LS,R) , AA(LS) ,BB(LS) ,FACT(7) 

DATA FACT/1 ,000 ,2,0DQ ,6 .000,2*1 ,0D0 , 

1 120. 0D0, 720, 0D0, 50*10. 0D0/ 


IF (TERM2 ,EQ, 0) GO TO 20 

COMPUTES S(INDEX+J ,I)sS(TERMl ,I)*“POWER » S(TERM2+J,I) 
FOR JsO,1,,..,TERM3-TERM2 

KK a ~1 

DO 10 K s TERM2,TERM3 

KK * KK + 1 

INDX s INDEX + KK 

AA(INDX) = AA(TERM1)**P0WER » AA(K) 

IPOW s POWER 

IF (TERM1 ,EQ, K) IPOW s IPOW + 1 

BB(INDX) a BB(TERM1)»»POWER * BB(K) “FACT(IPOW) 

DO 10 I s 2,R 

10 S(INDX,I) s S(TERM1,I)«F0WER * S(K,I) 

INDB*/ b INDEX + TERM3 - TERM2 + 1 


RETURN 

COMPUTES S( INDEX, I) b S(TERM1 ,I)“»POWER 
20 CONTINUE 

AA( INDEX) s AA(TERM1 ) ““POWER 
BB( INDEX) a BB(TERM1)“ “POWER * FACT(POWER) 
DO 25 I s 2, R 

25 S(INDEX, I) b S(TERM1 ,I)*“POWER 
INDEX = INDEX + 1 

C 

RETURN 

END 
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